home *** CD-ROM | disk | FTP | other *** search
/ MacWorld 1997 September / Macworld (1997-09).dmg / Serious Software / Cherwell Scientific Demos / pro Fit / pro Fit 5.0 demo (fpu).sea / pro Fit 5.0 demo (fpu) / Notes / Fitting differential equations / Example 1 / Nonlinear Decay function < prev    next >
Text File  |  1996-06-01  |  3KB  |  147 lines

  1. {
  2.  The following are the four integrating functions NonlinearDecay,
  3.  NonlinearDecay2, NonlinearDecay3, NonlinearDecay4, which are
  4.  described in the TechNote "Fitting differential equations".
  5. }
  6.  
  7.  
  8. function NonlinearDecay;
  9.  const step = 0.01;
  10.  var t, dt;
  11. begin
  12.  y := a[4];    {start value}
  13.  t := 0;       {start time}
  14.  dt := step;   {step width}
  15.  while t < x do
  16.  begin
  17.    if t+dt > x then dt := x-t;
  18.    y := y + dt*(y*(a[1] + y*(a[2] + y*a[3])));
  19.    t := t+dt;
  20.  end;
  21. end;
  22.  
  23.  
  24. function NonlinearDecay2; {Runge-Kutta}
  25.  const step = 0.05;
  26.  var t, dt, m1, m2, m3, m4, y1;
  27. begin
  28.  y := a[4];    {start value}
  29.  t := 0;       {start time}
  30.  dt := step;   {step width}
  31.  while t < x do
  32.  begin
  33.    if t+dt > x then dt := x-t;
  34.    m1 := y*(a[1] + y*(a[2] + y*a[3]));
  35.    y1 := y + dt*m1*0.5;
  36.    m2 := y1*(a[1] + y1*(a[2] + y1*a[3]));
  37.    y1 := y + dt*m2*0.5;
  38.    m3 := y1*(a[1] + y1*(a[2] + y1*a[3]));
  39.    y1 := y + dt*m2;
  40.    m4 := y1*(a[1] + y1*(a[2] + y1*a[3]));
  41.    y := y + dt*(m1 + 2*m2 + 2*m3 + m4)/6;
  42.    t := t+dt;
  43.  end;
  44. end;
  45.  
  46.  
  47. function NonlinearDecay3; {temporary storage}
  48.  const step = 0.01;
  49.  var t, dt;
  50.  var lastX, lastY;
  51.  
  52. procedure First; {called when lastX, lastY are stale}
  53. begin
  54.  lastX := 0;
  55.  lastY := a[4];
  56. end;
  57.  
  58. begin
  59.  if x >= lastX then begin
  60.    y := lastY;    {start value}
  61.    t := lastX;    {start time}
  62.  end
  63.  else begin
  64.    y := a[4];     {start value}
  65.    t := 0;        {start time}
  66.  end;
  67.  dt := step;   {step width}
  68.  while t < x do
  69.  begin
  70.    if t+dt > x then dt := x-t;
  71.    y := y + dt*(y*(a[1] + y*(a[2] + y*a[3])));
  72.    t := t+dt;
  73.  end;
  74.  lastX := x;
  75.  lastY := y;
  76. end;
  77.  
  78.  
  79.  
  80. function NonlinearDecay4; {temporary storage / Derivatives}
  81.  const step = 0.01;
  82.        epsilon = 1e-8; {a small deviation for all params}
  83.  var lastX0, lastY0, lastX1, lastY1, lastX2, lastY2, 
  84.      lastX3, lastY3, lastX4, lastY4;
  85.  
  86.  
  87. function Integrate(lastX, lastY, xVal, a1, a2, a3, a4);
  88.  {runs one integration}
  89.  var yVal, t, dt;
  90. begin
  91.  if xVal >= lastX then begin
  92.    yVal := lastY;    {start value}
  93.    t := lastX;    {start time}
  94.  end
  95.  else begin
  96.    yVal := a4;       {start value}
  97.    t := 0;        {start time}
  98.  end;
  99.  dt := step;   {step width}
  100.  while t < xVal do
  101.  begin
  102.    if t+dt > xVal then dt := xVal-t;
  103.    yVal := yVal + dt*(yVal*(a1 + yVal*(a2 + yVal*a3)));
  104.    t := t+dt;
  105.  end;
  106.  Integrate := yVal;
  107. end;
  108.  
  109.  
  110. procedure First; {called when lastX, lastY are stale}
  111. begin
  112.  lastX0 := 0; lastY0 := a[4];
  113.  lastX1 := 0; lastY1 := a[4];
  114.  lastX2 := 0; lastY2 := a[4];
  115.  lastX3 := 0; lastY3 := a[4];
  116.  lastX4 := 0; lastY4 := a[4]+epsilon;
  117. end;
  118.  
  119. procedure Derivatives;
  120.   var yVal, yVal1;
  121. begin
  122.   yVal := Integrate(lastX0, lastY0, x, a[1], a[2], a[3], a[4]);
  123.   lastX0 := x; lastY0 := yVal;
  124.  
  125.   yVal1 := Integrate(lastX1, lastY1, x, a[1]+epsilon, a[2], a[3], a[4]);
  126.   lastX1 := x; lastY1 := yVal1;
  127.   dyda[1] := (yVal1-yVal)/epsilon;
  128.  
  129.   yVal1 := Integrate(lastX2, lastY2, x, a[1], a[2]+epsilon, a[3], a[4]);
  130.   lastX2 := x; lastY2 := yVal1;
  131.   dyda[2] := (yVal1-yVal)/epsilon;
  132.  
  133.   yVal1 := Integrate(lastX3, lastY3, x, a[1], a[2], a[3]+epsilon, a[4]);
  134.   lastX3 := x; lastY3 := yVal1;
  135.   dyda[3] := (yVal1-yVal)/epsilon;
  136.  
  137.   yVal1 := Integrate(lastX4, lastY4, x, a[1], a[2], a[3], a[4]+epsilon);
  138.   lastX4 := x; lastY4 := yVal1;
  139.   dyda[4] := (yVal1-yVal)/epsilon;
  140. end;
  141.  
  142. begin
  143.  y := Integrate(lastX0, lastY0, x, a[1], a[2], a[3], a[4]);
  144.  lastX0 := x;
  145.  lastY0 := y;
  146. end;
  147.